Nuage d’ions¶
On considère un nuage de \(N\) ions de masse \(m\) et charge \(q\). On note \(\vec{r_i}\) et \(\vec{v_i}\) la position et la vitesse du ième ion. Les ions sont dans un piège électrostatique. De plus, ils interagissent entre eux par la force de Coulomb.
La force électrostatique dérive d’un potentiel électrostatique :
La force de Coulomb s’écrira sous la forme:
Pour simplifier et avoir une meilleure représentation graphique, on suprimera totalement la dimension \(z\).
On utilisera des unités adimensionnées avec des constantes de l’ordre de 1. Par défaut on prendra : \(m=9\), \(q=1\), \(k_x=1\), \(k_y=1.3\), \(\kappa=1\).
L’état des \(N\) particules à un instant donné est représenté par 4 tableaux numpy de taille \(N\): r_x, r_y, v_x, v_y. Toutes les fonctions seront écrites avec ces variables. On ne regroupera les 4 tableaux que dans la fonction f(t, y) qui sera utilisée par solve_ivp.
Le tableau y est défini en rassemblant les 4 tableaux numpy. On définit alors les fonctions join et split (voir ci-dessous)
Les paramètres seront des constantes globales.
Écrire la fonction
force_piege(r_x, r_y)qui renvoie la force dérivant du potentiel (f_xetf_y).Écrire la fonction
force_coulomb(r_x, r_y)qui renvoie \(f_x\) et \(f_y\), la force de Coulomb. On ne cherchera pas à éviter les boucles for dans cette fonction.Écrire la fonction
f(t, y)qui définit la dynamique du problème.On condidère un nuage avec une distribution initiale donnée par une loi de Maxwell-Boltzmann (avec \(k_B=1\) et \(T=1\)) pour les particules sans interaction. Calculer jusqu’au temps \(T=20\) l’évolution de la position des particules.
Vérifier que l’énergie totale est conservée. On utilisera les fonctions ci-dessous.
Faire une animation (voir exemple de code ci dessous)
On rajoute une force de dissipation, selon l’axe \(x\) :
\[ \vec F_{i,x} = -\alpha v_{i, x} \]Simuler l’experience en prenant \(\alpha = 0.1\) et \(N=20\). Que se passe-t-il au temps long ?
Utiliser le decorateur
jitde numba pour accélérer le code. Faire une comparaison de vitesse.On considère un nuage contenant deux types de particules (masse \(m_1\) et \(m_2\)) et de même charge \(q=1\). Seuls les particules de type \(1\) sont soumises à la force dissipative. Simuler l’expérience en prenant les mêmes paramètres qu’à la question 7. On prendra par exemple \(m_1=9\), \(N_1=20\) et \(m_2=2\) et \(N_2=5\)
Solution¶
import matplotlib.pyplot as plt
import numpy as np
from numpy import pi
from scipy.integrate import solve_ivp
# Paramètres
m = 9; q = 1
k_x, k_y = 1, 1.3
kappa=1
k_B = 1
T = 1
Quelques fonctions pour vous aider :
Fonction pour passer de 4 tableaux à un seul et réciproquement
def join(r_x, r_y, v_x, v_y):
return np.concatenate((r_x, r_y, v_x, v_y))
def split(y):
N = len(y)//4
return y[:N], y[N:2*N], y[2*N:3*N], y[3*N:4*N]
Calcul de l’énergie
def energie_cinetique(r_x, r_y, v_x, v_y):
return np.sum(.5*m*v_x**2 + .5*m*v_y**2)
def energie_piege(r_x, r_y):
return np.sum(k_x*r_x**2/2 + k_y*r_y**2/2)
def energie_coulomb(r_x, r_y):
N = len(r_x)
total = 0
for i in range(N-1):
for j in range(i+1, N):
d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2
total += kappa*q**2/np.sqrt(d2)
return total
def energie_totale(r_x, r_y, v_x, v_y):
return (energie_coulomb(r_x, r_y) +
energie_piege(r_x, r_y) +
energie_cinetique(r_x, r_y, v_x, v_y))
# Question 1
def force_piege(r_x, r_y):
return -k_x*r_x, -k_y*r_y
# Question 2
def force_coulomb(r_x, r_y):
N = len(r_x)
f_x = np.zeros(N)
f_y = np.zeros(N)
for i in range(N):
tot_x, tot_y = 0, 0
for j in range(N):
if i!=j:
d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2
coef = kappa*q**2/(np.sqrt(d2)**3)
tot_x += coef*(r_x[i]-r_x[j])
tot_y += coef*(r_y[i]-r_y[j])
f_x[i] = tot_x
f_y[i] = tot_y
return f_x, f_y
# Question 3
def force_totale(r_x, r_y, v_x, v_y):
piege = force_piege(r_x, r_y)
coulomb = force_coulomb(r_x, r_y)
return piege[0]+coulomb[0], piege[1]+coulomb[1]
def f(t, y):
r_x, r_y, v_x, v_y = split(y)
f_x, f_y = force_totale(r_x, r_y, v_x, v_y)
return join(v_x, v_y, f_x/m, f_y/m)
# Question 4
N = 50
r_x_0 = np.random.normal(size=N)
r_y_0 = np.random.normal(size=N)
v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
t_max = 20
t_eval = np.linspace(0, t_max, 51)
res = solve_ivp(f, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0),
t_eval=t_eval)
# Question 5
energie = []
for i, t in enumerate(t_eval):
r_x, r_y, v_x, v_y = split(res.y[:,i])
energie.append(energie_totale(r_x, r_y, v_x, v_y))
energie = np.array(energie)
plt.plot(t_eval, energie)
[<matplotlib.lines.Line2D at 0x7f10b775e8b0>]
* code pour afficher une animation. res est le tableau provenant de solve_ivp
from IPython.display import HTML
from matplotlib.animation import FuncAnimation
fig, ax = plt.subplots()
l, = ax.plot(split(res.y[:,0])[0], split(res.y[:,0])[1], 'o')
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
def animate(i):
r_x, r_y, v_x, v_y = split(res.y[:,i])
l.set_data(r_x, r_y)
ani = FuncAnimation(fig, animate, frames=len(res.t))
HTML(ani.to_jshtml())
#Question 7
alpha = 0.2
def force_friction(r_x, r_y, v_x, v_y):
return -alpha*v_x, 0*v_y
def force_totale(r_x, r_y, v_x, v_y):
piege = force_piege(r_x, r_y)
coulomb = force_coulomb(r_x, r_y)
friction = force_friction(r_x, r_y, v_x, v_y)
return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]
def f(t, y):
r_x, r_y, v_x, v_y = split(y)
f_x, f_y = force_totale(r_x, r_y, v_x, v_y)
return join(v_x, v_y, f_x/m, f_y/m)
N = 20
r_x_0 = np.random.normal(size=N)
r_y_0 = np.random.normal(size=N)
v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
t_max = 1500
t_eval=np.linspace(0, t_max, 400)
%%time
res = solve_ivp(f, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0),
t_eval=t_eval)
CPU times: user 10 s, sys: 59 µs, total: 10 s
Wall time: 10 s
energie = []
for i, t in enumerate(t_eval):
r_x, r_y, v_x, v_y = split(res.y[:,i])
energie.append(energie_cinetique(r_x, r_y, v_x, v_y))
energie = np.array(energie)
plt.semilogy(t_eval, energie)
[<matplotlib.lines.Line2D at 0x7f10b450fbb0>]
fig, ax = plt.subplots()
l, = ax.plot(split(res.y[:,0])[0], split(res.y[:,0])[1], 'o')
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
def animate(i):
r_x, r_y, v_x, v_y = split(res.y[:,i])
l.set_data(r_x, r_y)
ani = FuncAnimation(fig,
animate,
frames=len(res.t),
interval=100)
HTML(ani.to_jshtml());
# Question 8
from numba import jit
alpha = 0.2
@jit
def join_numba(r_x, r_y, v_x, v_y):
return np.concatenate((r_x, r_y, v_x, v_y))
@jit
def split_numba(y):
N = len(y)//4
return y[:N], y[N:2*N], y[2*N:3*N], y[3*N:4*N]
@jit
def force_piege_numba(r_x, r_y):
return -k_x*r_x, -k_y*r_y
@jit
def force_coulomb_numba(r_x, r_y):
N = len(r_x)
f_x = np.zeros(N)
f_y = np.zeros(N)
for i in range(N):
tot_x, tot_y = 0, 0
for j in range(N):
if i!=j:
d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2
coef = kappa*q**2/(np.sqrt(d2)**3)
tot_x += coef*(r_x[i]-r_x[j])
tot_y += coef*(r_y[i]-r_y[j])
f_x[i] = tot_x
f_y[i] = tot_y
return f_x, f_y
@jit
def force_friction_numba(r_x, r_y, v_x, v_y):
return -alpha*v_x, 0*v_y
@jit
def force_totale_numba(r_x, r_y, v_x, v_y):
piege = force_piege_numba(r_x, r_y)
coulomb = force_coulomb_numba(r_x, r_y)
friction = force_friction_numba(r_x, r_y, v_x, v_y)
return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]
@jit
def f_numba(t, y):
r_x, r_y, v_x, v_y = split_numba(y)
f_x, f_y = force_totale_numba(r_x, r_y, v_x, v_y)
return join_numba(v_x, v_y, f_x/m, f_y/m)
N = 20
r_x_0 = np.random.normal(size=N)
r_y_0 = np.random.normal(size=N)
v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)
t_max = 1500
t_eval=np.linspace(0, t_max, 400)
# Ici, on apelle f_numba un première fois pour la compilation
f_numba(0.0, join(r_x_0, r_y_0, v_x_0, v_y_0));
%%time
res = solve_ivp(f_numba, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0),
t_eval=t_eval)
CPU times: user 168 ms, sys: 0 ns, total: 168 ms
Wall time: 167 ms
# Question 9
particle_type = np.array([0]*20 + [1]*5)
masses = np.array([9, 3]) # The masse of particles given by masses[particle_type]
colors = np.array(['blue', 'red'])
@jit
def force_friction_numba(r_x, r_y, v_x, v_y, particle_type):
return -alpha*v_x*(particle_type==0), 0*v_y
@jit
def force_totale_numba(r_x, r_y, v_x, v_y, particle_type):
piege = force_piege_numba(r_x, r_y)
coulomb = force_coulomb_numba(r_x, r_y)
friction = force_friction_numba(r_x, r_y, v_x, v_y, particle_type)
return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]
@jit
def f_numba(t, y, particle_type):
r_x, r_y, v_x, v_y = split_numba(y)
f_x, f_y = force_totale_numba(r_x, r_y, v_x, v_y, particle_type)
return join_numba(v_x, v_y, f_x/masses[particle_type], f_y/masses[particle_type])
N = len(masses[particle_type])
r_x_0 = np.random.normal(size=N)
r_y_0 = np.random.normal(size=N)
v_x_0 = np.random.normal(scale=1/np.sqrt(masses[particle_type]), size=N)
v_y_0 = np.random.normal(scale=1/np.sqrt(masses[particle_type]), size=N)
t_max = 1500
t_eval=np.linspace(0, t_max, 400)
res = solve_ivp(f_numba, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0),
t_eval=t_eval, args=(particle_type,))
fig, ax = plt.subplots()
X = split(res.y[:,0])[0]
Y = split(res.y[:,0])[1]
l = ax.scatter(X, Y, c=colors[particle_type])
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
def animate(i):
r_x, r_y, v_x, v_y = split(res.y[:,i])
off = np.array([r_x, r_y]).T
l.set_offsets(off)
ani = FuncAnimation(fig,
animate,
frames=len(res.t),
interval=100)
HTML(ani.to_jshtml())